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I. INTRODUCTION 

Random systems represent the paradigm of spatially 
inhomogeneous systems where scale invariance is only 
preserved on average, but not for specific disorder re- 
alizations [0. In such systems, not a single exponent but 
instead an infinite hierarchy of independent exponents 
are expected to describe the scaling behaviour of local 
quantities like order parameter density profiles or corre- 
lation functions. This property, linked to the non self- 
averaging behaviour of the corresponding physical quan- 
tities is ususally referred to as multifractality. The 
keystone concept which enters the description of multi- 
scaling properties is that of scaling dimensions associated 
to the moments of the local physical property, or equiva- 
lently the universal function H(a) corresponding to the 
Legendre transform of the set of independent scaling in- 
dexes (See e.g. Refs. ||J^). 

Ten years ago already, in a series of illuminating pa- 
pers, Ludwig |8|J^ and Ludwig and Cardy [p^ reported 
an extensive analytic study of 2D random bond Potts 
ferromagnets in the regime where bond randomness is 
slightly relevant, q close to 2, g being the number of states 
per spin. Their studies included perturbative calcula- 
tions of the conformal anomaly, of the thermal scaling 
dimension and of the multifractal behaviour of spin-spin 
correlation functions. 

The essential of the numerical studies dealing with 
scaling dimensions of average quantities in random bond 
Potts models at small values of q were performed by 
Monte Carlo (MC) simulations combined to standard Fi- 



nite Size Scaling (FSS) or transfer matrix (TM) 

calculations associated to conformal methods |p5|-p^ at 
the random fixed point (FP) of self-dual disordered mod- 
els. In particular, an excellent quantitative agreement for 
the magnetic scaling dimension in the three-state Potts 
model was reported in Refs. In which concerns 

the multiscaling properties of spin-spin correlation func- 
tions, although some of Ludwig's predictions have par- 
tially been verified both in cylinder geometry [|l5|,^ and 
in square geometry ||2l|] , the agreement with analytical 
expansions was less conclusive and in particular the shape 
of the probability distribution has not been reproduced. 

Monte Carlo simulations are not convenient to study 
numerically the vicinity of the Ising model, q — 2, where 
perturbation expansions are supposed to apply. The 
number of states per spin, q, is indeed restricted to inte- 
ger values in such simulations. In this paper, we therefore 
use a rather different approach already used by different 
authors |l5yi^ , [l9| , ^ . This technique is based on the 
Fortuin-Kasteleyn graph representation which en- 
ables TM calculations (24| where q enters as a parameter 
that can take non integer values. We also benefit from 
previous studies with a bimodal probability distribution 
of spin-spin interactions where the disorder amplitude 
was found to have a deep influence on the measured crit- 
ical properties in numerical studies It should thus 
be chosen carefully in order to avoid crossover perturba- 
tions due to the unstable pure model fixed point, r = 1, 
and the percolation fixed point, r — > 00, r being the ratio 
between strong and weak couplings. 

In this paper, we are mainly interested in the multi- 
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scaling behaviour of the spin-spin correlation functions. 
The first section reminds the reader of the essential rel- 
evant theoretical results which have been obtained by 
several groups using perturbative techniques around the 
pure models conformal field theories. Section ^ explains 
the methodology and section IV gives the numerical re- 
sults: 

i) Universality in quenched disordered ferromagnetic 

Potts models is checked using different types of 
probability distributions of nearest-neighbour cou- 
plings. 

ii) The numerical study of the decay exponent of the sec- 

ond moment of the spin-spin correlation function 
is then compared to perturbative results in order 
to test Replica Symmetry and Replica Symmetry 
Breaking scenarios. 

iii) Finally, we study other moments and deduce the 

shape of the universal functions H{a) for different 
q values. The probability distribution of spin-spin 
correlation functions is also analyzed. 

The values of q are chosen in the range 2 to 4, where the 
pure model exhibits a second order phase transition, with 
a special attention paid on the neighbourhood of q = 3. 



II. SUMMARY OF THE PERTURBATION 
RESULTS 

A. The 2D random Ising model 



denotes the disorder average. The central charge c, de- 
fined by the leading size dependence of the free energy 
density in the cylinder geometry, ^% , thus exhibits loga- 
rithmic corrections which make its exact value i difficult 
to extract numerically 

Using a perturbation expansion, Ludwig later obtained 
the behaviour of the moments of the spin-spin correlation 
function [01: 



(2) 



when A, the strength of disorder, is strong enough. 
Brackets denote the thermal average. We can also intro- 
duce a reduced correlation function whose leading power- 
law behaviour is 



(3) 



Ludwig's results imply that logarithmic corrections are 
absent in the case of the average correlation function (p — 
1), {a{0)a{p)) ^ P~^^^, while typical correlation func- 
tions (p = 0) exhibit such corrections, exp hi{a{0)a{p)) ^ 
p-i/4^^}j^p~)-i/8^ observed numerically Eo]. Further- 
more, a unique scaling dimension Xa — describes the 
leading power-law decay of all the moments of the spin- 
spin correlation function, {a{0)a{p))P ^ p''^^^" , and no 
multiscaling behaviour is expected apart from the loga- 
rithmic correction term. 

The surface correlation function has also been studied 
recently and was found to be self- averaging p| . 



According to the celebrated Harris criterion p5| , 
quenched randomness is a marginal perturbation in the 
2D Ising model. This situation has focused a consid- 
erable interest on the critical properties of the random 
bond Ising model (RBIM), and, after partially conflicting 
results, disorder was eventually found to be marginally 
irrelevant, leading to an unchanged universal behaviour 
apart from logarithmic corrections for the ensemble av- 
erage of some physical quantities p^-|3l|]. These results 
were then carefully checked through intensive MC simula- 
tions p^-|3^ and series expansions ps]-^. We mention 
that site dilution is still subject to controversial inter- 
pretations (See e.g. Ref. [^) in spite of a conclusive 
recent work leading to the same conclusions than bond 
randomness p^ . 

The free energy density in a strip of width L with pe- 
riodic boundary conditions was obtained by Ludwig and 
Cardy 0: 
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(1) 



where A, the variance of the Gaussian probability dis- 
tribution of exchange interactions, is the strength of dis- 
order related to the ratio r. In this expression, overbar 



B. The random bond Potts model 

In the case of the Potts model with q > 2, disorder 
is a relevant perturbation which modifies the universal 
critical behaviour and leads to new fixed point critical 
properties. Studying the effect of a slightly relevant per- 
turbation on the finite-size scaling behaviour of the free 
energy density in a strip geometry of width L at the new 
fixed point, Ludwig and Cardy |L0||, obtained perturba- 
tively the central charge c'{q) of the 2D g— state Potts 
model with weak quenched bond randomness (here and 
in the following, primes denote the central charge and the 
critical exponents at the disordered fixed point, while un- 
primed symbols refer to the pure fixed point quantities). 
Using an expansion in g — 2 around the Ising model, the 
random anomaly was deduced from the random free en- 
ergy of the strip. 



fiL,g*)^A-^ + 0{L-% 



(4) 



given in the replica formalism by the quenched free en- 
ergy df{n)/dn\.n^o- 



c'{q) 
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where yn is the renormaUzation group (RG) eigenvalue 
associated to the bond disorder |2^, yH = a/v = 
2 — 2X5(5), ^eil) being the scahng dimension of the en- 
ergy density in the pure model ^. This latter dimen- 
sion is obtained for arbitrary g < 4 by the den Nijs con- 
jectureE^,^, rigourously proved by Dotsenko and Fa- 
teev ill: 



1 



< ^ = — cos 

TT 



(iV^) < 1, (6) 



and to lowest order, the RG eigenvalue is proportional to 
q-2: yH = ^{q-1) + 0[{q-2f]. The deviation of the 
conformal anomaly from its pure fixed point value, 



c{q) = 1 



A* 



(7) 



is difficult to measure, since it is only of third order in 
yn and it requires a very good accuracy to distinguish 
between pure and random values. 

The thermal exponent was similarly obtained to two- 
loop order by Ludwig 
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x'M = xM + -yH + ^Vh + 0{y]i) 



= 1 



Oiy'n), 



(8) 



and the three-loop correction was reported by Jug and 
Shalaev [||. 

The correction to the magnetic exponent requires a 
three-loop calculation. It was obtained by Dotsenko et 
al |6| 

xUq) = x^q) + ^^^T^^/h + 0{y%), (9) 

and checked by Picco using MC simulations . In con- 
tradistinction with the thermal exponent, the deviation 



from the pure fixed point value is quite small close to 
g = 2. 

Searching for multiscaling properties, Ludwig ob- 
tained, up to linear order, the scaling dimension of the 
p*''-moment ^ p| of the reduced spin-spin correlation 



-i/p 



22:„p(<?) g^j-^j Lewis per- 



function, {a{Q)a{p))P ~ p 
formed recently the computation up to the second or- 
der MM 



B{p~2)]y% + 0{yj,), (10) 



where A = ii - 41n2 and B = 5^(33 _ 29\/37r/3). Here, 
the exponent corresponding to the average critical corre- 
lation function at the random fixed point is denoted by 
x' i{q) = x'^{q), while the typical behaviour corresponds 
to p = 0. This result, obtained in the Replica Symme- 
try (RS) scenario, contains the special case of the second 
moment performed by Dotsenko et al in order to 
compare between Replica Symmetry 



x'^2iq) = x^{q) 
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16 



yn 



If41n2-ii 
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Vh + Oiy'n), 



(11) 



and Replica Symmetry Breaking (RSB) 



x'a^iq) 



Xa{q) - Y^yH 
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^In this paper, we use the notations of Refs. p|-p^: The 
expansion parameter is i/h- It is related to the parameter 
e, linked to the deviation from the pure Ising model central 
charge, used by Dotsenko and co-workers in Refs. p^ , p7|j4£| , |50[ ] 
by yn = 3e. 

^In the literature, many different notations have been used 
for the scaling dimensions of the moments of the correlation 
functions. Our notation corresponds to that of Lewis | ^ : 
x'^p <-> A'^p. The correspondence with other works is the 
following: Ludwig x'^p « Xn/N, Dotsenko et al |p9| : 
x'^2 ^ A^2/2 and Olson and Young |2lj: x'^p <-» ri„/2. 

^InRef. 0, the thermal and magnetic exponents have been 
computed with both RS and RSB scenarios. While Eq. (^) 
for the average behaviour is unchanged up to the third order 
in the RSB scheme, the thermal exponent in Eq. (^) becomes 
x':ig) = l + 0{yj,). 
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III. MODEL AND METHODOLGY 



In this paper, we consider Potts-spin variables, <7j e 
1,2, ... ,q on the sites of a square lattice with indepen- 
dent quenched random nearest-neighbour ferromagnetic 
interactions Kij. These exchange couplings are taken 
from a probability distribution V{Kij), and in most of 
our applications, they can take two values, Ki = K, 
K2 — rK, r > 1 with equal probabilities: 



1 



1 



ViK,,) = -S{K,, -K) + -5{K,, - rK), 



(13) 



where r measures the strength of disorder. Our method- 
ology will be discussed in this section with the particular 
case of probability distribution (|l^) and the generaliza- 
tion to other distributions will be presented in the next 
section. 

The Hamiltonian of the model is thus written 



(14) 



We only consider self-dual models for which the crit- 
ical point is exactly known. In the case of the bimodal 
distribution (p^), the self-duality point 



[exp(i^,(r)) - l][exp(ri^,(r)) - I] = q 



(15) 



corresponds to the critical point of the model if only one 
phase transition takes place in the system as rigourously 
shown in Ref. [ pT| . 

The degree of dilution in the system can be varied by 
changing the ratio of the strong and weak couplings, r. 
At r = 1, one recovers the perfect q-state Potts model, 
whereas for r ^ cx) we are in the percolation limit, where 
Tc = 0. The intermediate regime of dilution \ < r < 00 
is expected to be controlled by the random fixed point 
located at some r — r*{q) This optimal disorder 

amplitude, r*{q), can be obtained numerically from the 
maximum condition of the effective central charge of the 
disordered system. Dotsenko and co-workers for example 
considered n g— state Potts models coupled via energy- 
energy interactions and obtained perturbatively the cen- 
tral charge deviation from the decoupling limit (where c 
is given by the sum of the central charges of the decoupled 
models) ||: Ac = -^jj^y^ + 0{y%). For n > 1, 
Ac satisfies the Zamolodchikov's c— theorem according 
to which there exists a c— function decreasing along RG 
flows and giving the central charge at the fixed point . 
In the case of random systems (n —t in the replica 
approach), the central charge increases and can be ex- 
pected to reach a maximum value at an optimal disorder 
amplitude where the random FP exponents may be ex- 
tracted from numerical data. This property, linked to 
non-unitarity in the presence of disorder, is indeed ob- 
served in simulations |l^Jl^,^ . 

In the following we used a TM technique, based on the 
Blote and Nightingale connectivity transfer matrix pi]. 



which enables to compute the physical quantities in long 
cylinders. Since transfer operators in the time direction 
do not commute in disordered systems, the free energy 
density is defined by the leading Lyapunov exponent. For 
an infinitely long strip of width L with periodic bound- 
ary conditions, the leading Lyapunov exponent is given 
by the Furstenberg method p% : 



Ao(L) = lim — In 

m^oo rn 



\k=l 



(16) 



where is the transfer matrix and | wo) is a unit initial 
vector. The quenched free energy density is thus given 
by 



/(L) = -L-iAo(L). 



(17) 



The shape of the central charge as a function of the 
disorder amplitude is shown in Fig. Q for several values 
of q. Each realization of disorder is obtained via 10^ iter- 
ations of the transfer matrix and the free energy density 
was averaged over 96 such realizations. 



^=1.5 



q=2.5 ^=3 



1 \ ' 



1^ 




0.9998 



FIG. 1. Behaviour of the central charge as a function of 
the disorder amplitude for different values of the number of 
states q (binary distribution of Eq. (jl^)). The solid lines 
are parabolic fits. The maximum corresponds to the optimal 
value of disorder amplitude. 
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TABLE I. Optimal disorder amplitude and corresponding 
values of the central charges of the disordered Potts model for 
different values of q for the binary probability distribution in 
Eq. (p^. The last columns give the numerical values of the 
pure model central charge c(g) (Eq. (^) and of the expan- 
sion parameter used in the perturbation resuhs, and shows 
that the perturbation expansion obviously breaks down when 
q increases. 



q 


r*{q) 


TM resuh 


Eq. i) 


c{q) 




1.5 


2.25 


0.283(1) 




0.288 




2. 


2.49 


0.496(6) 


0.5 


0.500 


0. 


2.25 


3.18 


0.584(8) 


0.5876 


0.588 


0.1036 


2.5 


3.96 


0.662(8) 


0.6661 


0.666 


0.2036 


2.75 


4.70 


0.732(8) 


0.7371 


0.736 


0.3017 


3. 


5.36 


0.797(8) 


0.8020 


0.800 


0.4000 


3.25 


5.94 


0.857(9) 


0.8617 


0.858 


0.5013 


3.5 


6.45 


0.912(11) 


0.9174 


0.910 


0.6101 


4. 


7.30 


1.011(12) 


1.0312 


1.000 


1. 



The central charge follows from a polynomial fit 



I{L) 



AaL" 



AaL- 



AhL- 



(18) 



where the remaining coefficients Ai have been included 
in order to simulate finite size effects, although the exact 
dependence of corrections to scaling is not known. The 
calculations were performed on strips of widths L = 2 to 
8. The central charge is still sensible to the number of 
disorder configurations entering the average and to the 
degree of the polynomial fit. Nevertheless, it is expected 
that the position of the maximum presents small enough 
deviations for the critical exponents to reflect the dis- 
ordered fixed point regime in the neighbourhood of this 
maximum. We have checked different types of fits, with 
less parameters and also including logarithmic terms in 
the vicinity oi q — 2 where equation (m) is supposed to be 
valid, but we were not able to improve the results, since, 
as it has already been observed the values of c'(g) are 
systematically below the perturbative result of Ludwig 
and Cardy ||l^, and even below the pure model central 
charge at small values of q. In all the cases, the error can 
be estimated by the fluctuations of the results with dif- 
ferent fitting procedures and it is of order 10~^ to 10"^. 
For example we obtain c'(2) = 0.496 with Eq. (|l|), 0.492 
with Ad =^ As ^ 0, 0.497 with As = but a (Ini)-^ 
term added or 0.495 with Aq = As = and the log-term 
present. The error bars given in table | correspond to 
the fluctuations of the results obtained with the different 
fits. The maximal values and optimal disorder amplitude 
r* (q) are also given in table |. 

For a specific disorder realization, the spin-spin corre- 
lation function along the long direction, w, of the strip 



(19) 



where (. . .) denotes the thermal average, is given by a 
product of non-commuting transfer matrices. They were 
computed on strips of widths L = 2 to 8 and then aver- 
aged over 80 000 disorder realizations. 

We will now assume that conformal covariance can be 
applied to the order parameter correlation function and 
its moments. In the infinite complex plane z = x + iy the 
correlation function and its moments exhibit the usual 
algebraic decay at the critical point 



{a{zi)a{z2))p'^'' - const x p-2<p('?), 



(20) 



where p =\ zi — Z2\. Multiscaling arises when the expo- 
nents x'^p{q) are all different, depending upon the value 
of the moment order p, and their p— dependence is a con- 
vex function ||^. Under the logarithmic transformation 
w = -^Inz = u + iv which maps the infinite plane geom- 
etry inside an infinitely long strip of width L, one gets 
the exponential decay along the strip 



rl/p 



{a{j)a{j + u))P = const x cxp 



2tt 



(21) 



The scaling dimensions x'^p (q) at different strip sizes can 
thus be deduced from an exponential fit, and a quadratic 
extrapolation at i ^ oo is performed to get the corre- 
sponding value in the thermodynamic limit. The calcu- 
lation of errors follows the lines explained in Rcf. |l9| ]. 

This method was used in Ref. ||l^ for the average cor- 
relation function, i.e. for p = 1, and in the large q regime. 
In this work, we calculate the higher moments, as well as 
the typical behaviour, governed by the derivative of the 
exponent x'^p with respect to p, evaluated in the limit 
p — > 0. In previous works, it was shown that the re- 
sults from exponent extrapolations are sensitive to the 
value of the disorder strength chosen for the simulation, 
since the finite size corrections are very strong unless 
the calculations are performed close to the random fixed 
point as obtained from the maximum condition of 
the central charge of the model [ p^j52| . Our simulations 
were performed at these fixed point values of the disor- 
der strength, but for comparison we have also considered 
systems with somewhat different values. 

Since the correlation functions are not self-averaging, 
the disorder average must be performed carefully. We 
follow the same procedure as in Ref. where we com- 
pared the ensemble average {a{0)a{u)) to a cumulant ex- 
pansion in terms of the moments of ln((7(0)(T(u)), which 
are self- averaging M: 



(cr(0)cr(w))P = exp pln(cr(0)cr(u)) 

+ ip2 (in(a(0)a(u))2-ln(a(0)a(u))') + ... 

Although the average should in principle be done using 
the cumulant expansion, we observe that the direct av- 
erage, which is compatible with the cumulant expansion, 



(22) 
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is more stable than this latter expansion. This is partic- 
ularly true at high moment orders and is probably due 
to the large number of disorder realizations used in the 
calculations. This is illustrated for several moments at 
g = 3 in Fig. || where the solid lines represent the direct 
average over 96 000 different samples, while the open 
symbols correspond to the cumulant expansion up the 
the fifth order. 
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FIG. 2. Moments of the spin-spin correlation function 
p = 1, 2 and 3 from upper to lower curves {q = 3, L = 7, 
r = 5.36). The solid lines correpond to the average over 96 
000 different disorder realizations and the symbols are de- 
duced from the cumulant expansion up to the fifth order. 
The fluctuations become extremely large above the fifth or- 
der. Both solid lines and symbols give the same order for 
the corresponding scaling dimensions (related to the slopes 
of these curves), but the direct average leads to more precise 
results. 
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FIG. 3. Scaling dimension of the order parameter (binary 
disorder) compared to the third order expansion of Dotsenko 
and co-workers The scaling dimension corresponding to 
the pure model is shown for comparison. 

From the exponential decay of the average correlation 
function, the exponent x'^{q) is deduced and presented 
as a function of q for the case of a binary disorder in 
Fig. ^. Although these results are not new, since the 
same type of curve was reported by Cardy and Jacob- 
sen in Ref. |l8[ , the agreement with our results confirms 
the reliability of the averaging procedure. The results 
are compared to the third order expansion of Dotsenko 
et al [k 
table I 



6[ in Eq. (^). The data themselves are given in 
r The agreement is extremely good especially in 
the region where the expansion is supposed to be valid 
when q is not too far from the Ising model value q = 2. 

TABLE II. Comparison of the numerical results for the 
magnetic scaling dimension (bimodal probability distribu- 
tion) x'^{q) with the third order expansion of Dotsenko and 
co-workers The error bars systematically contain the 

analytical value. 



q 




T 


Expansion (B) 


TM result 


2 


0.12500 


0.1252(3) 


2.25 


0.12800 


0.1282(5) 


2.5 


0.13051 


0.1307(7) 


2.75 


0.13269 


0.1328(9) 


3. 


0.13465 


0.1347(11) 


3.25 


0.13653 


0.1364(13) 


3.5 


0.13845 


0.1379(14) 



Even close to the marginally irrelevant case of the Ising 
model where logarithmic corrections are known to be 
present for some quantities, we note that the numeri- 
cal data are quite satisfactorily in agreement with the 
perturbative results. This is due to the absence of loga- 
rithmic corrections for the average correlation function at 
q — 2, and we will see that this observation is no longer 
true in the following study of other moments. 



IV. TESTS OF UNIVERSALITY 

A. Universality of the average behaviour 

The question of universality in random systems is not 
yet solved, especially when frustration occurs, like in ran- 
dom fields or spin glasses [|5| . We address in this section 
the question of the influence of the particular shape of the 
probability distribution of exchange couplings on the uni- 
versality class. There are still subsisting doubts concern- 
ing universality since incompatible estimations (taking 
error bars into account) of the critical exponent x'^ were 
obtained with different probability distributions. For ex- 
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ample at g = 8, the choice of a bimodal probabihty distri- 
bution led to 0.151(4) IHl (FSS) or 0.1505(3) |9| (con- 
formal invariance), while a continuous distribution gave 
0.161(3) [HI (FSS). 

In this section, we show that the discrepancy is sim- 
ply due to crossover effects but does not imply absence 
of universality. Following the methodology of the previ- 
ous sections, the disordered fixed point regime is located 
at the maximum of the central charge and the scaling 
dimension x'^ of the average spin-spin correlation func- 
tions is estimated. The results of Fig. |^ are considered as 
a reference and the same method is applied to ternary, 
quaternary and continuous distributions and is shown to 
lead to identical critical exponents, within error bars, to 
those of a binary distribution. 



1. Ternary distribution 



The ternary probability distribution is defined by 
ViK.j) = i [6{K,j - Ko) + 5{K,j - K) + 5{K,, - rK)] 



with the self-duality condition 

[exp(K,(r)) - l][exp(rif,(r)) - 1] - [exp(i^o) - 1]' 



(23) 

= 9, 
(24) 

where Kq ~ ln(l + y/^) is the critical coupling of the 
pure system. The homogeneous system corresponds to 
the value r — 1. We tried different kinds of interpolation 
procedures for the free energy, as in Eq. (p^, but always 
found a monotonic variation of the corresponding central 
charge with respect to the disorder amplitude r. The ab- 
sence of maximum might be the sign of the presence of 
strong corrections, possibly due to the fact that one third 
of the exchange couplings keeps their pure value Kq even 
in the infinite-disorder limit. Nevertheless, we present in 
the table III the magnetic exponent as extracted from 
the average spin-spin correlation functions for different 
values of r. 

TABLE III. Scaling dimension x'^ of the average order pa- 
rameter for the q — 3 Potts model with a ternary distribution 
compared to the results for a binary distribution at the opti- 
mal disorder amplitude and for the pure model. 



Distribution 


disorder amplitude 


<(3) 


Pure 


r = 1 


0.1333 


Binary 


r* ~ 5.363 


0.1347(11) 


Ternary 


r = 2 


0.1339(3) 




r = 4 


0.1343(6) 




r = 5.363 


0.1344(8) 




r = 8 


0.1344(10) 




r = 12 


0.1343(13) 




r = 20 


0.1341(15) 



For strong disorder, the scaling dimension x'^, as pre- 
sented in the table III , shows a plateau with a value com- 
patible within error bars with that of the binary distri- 
bution, but the agreement is not yet conclusive, since the 
effective central charge was not found to display a clear 
maximum. 



2. Quaternary distribution 

The quaternary probability distribution is defined by 

r{K,j) = i [6{K,j -K) + 5{K,, - rK) 

+5{K,, - K') + 5{K,, - r'^K')] (25) 

where the four equi-probable exchange couplings if, rif , 
K' and r^K' are related by 

[exp(if,(r)) - l][exp(rif,(r)) - 1] 

= [cMKir)) - l]bMr^K{r)) - 1] = 9 (26) 

at the self-dual point of the model. The value r — 1 cor- 
responds to the pure system and the limit r +oo to a 
percolative regime. 

As seen on figure Fig. ^, the central charge presents 
both a maximum at r* ~ 2.000 and a "minimum" at 

~ 3.763. According to Zamolodchikov's c-theorem in 
the case of non-unitary theories, the latter case should 
correspond to an instable fixed point while the former 
situation is likely to be the disordered fixed point. 



In table IV, we collect the scaling dimensions of the 
average spin-spin correlation functions at these two fixed 
points and, for comparison, those of the pure model and 
of the disordered system with a binary distribution. 




FIG. 4. Behaviour of the central charge as a function of the 
disorder amplitude r for the quaternary distribution (^sj) at 
q = 3. The dotted curve is a guide for the eyes. 
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TABLE IV. Scaling dimensions x'„ {q — 3) of the average 
order parameter for the two extrema (r* (written in bold face) 
and r*) of the central charge with a quaternary distribution 
compared to the results for the binary case and for the pure 
model. 



Kij = In 



Distribution 


disorder amplitude 




Pure 


r = 1 


0.1333 


Binary 


r* ~ 5.363 


0.1347(11) 


Quaternary 


r* ~ 2.000 


0.1343(6) 




r, = 3.763 


0.0097(12) 



Inspection of the results in table IV reveals that com- 
patible exponents are obtained at the maximum of the 
central charge with both probability distributions. These 
data also confirm once a gain that the disorder amplitude 
r plays an essential role since the value at is def- 
initely excluded by the perturbative result. 



3. Continuous distribution 

Recently, Olson and Young |^ proposed slightly dif- 
ferent numerical estimations of the critical exponents of 
the average magnetization and of its first moments, com- 
pared to other independent studies [po| , ^^|j48| ) . They 
used an interesting continuous probability distribution 
of exchange couplings that they claimed to be less sen- 
sible to crossover effects. We show in the following that 
the randomness amplitude chosen in the simulations of 
Olson and Young is not the optimal one, since it is not 
strong enough to reach the disordered fixed point. This 
observation may be at the origin of the slight discrepancy 
between the extrapolated values of the exponents. 

The continuous probability distribution used by Olson 
and Young is generalized by introduction of a parame- 
ter r — which controls the strength of randomness. 
Following their notation, we have chosen the distribution 



where 



ny^. 



cosh ^ 



,K., _ 1 



(27) 



(28) 



Self-duality is ensured by the parity of the distribution 
ViVij)- The definition of the disorder amplitude r is such 
that r = 1 corresponds again to the pure system and 
r = e to the Olson- Young distribution. The probability 
distribution of exchange couplings Kij is given by 

r{K,j) = (29) 
and can be generated by the formula 



tan^ 



(30) 



if a; G [0; 1[ is a uniformly distributed random variable. 
Examples of probability distributions at different disor- 
der amplitudes are shown in Fig. |^. 



1.5 



r=1.5 

/ \ r=e (Olson- Young) 

/ \ r*=3.88I 

' r=6 



0.5 




FIG. 5. Probability distribution of exchange couplings 
V{Kij) for disorder amplitudes r in the range [1.5; 6]. 



0.79 




0.788 



0.786 - 



0.784 



0.782 



FIG. 6. Central charge of the g = 3 random bond Potts 
model for the continuous distribution Eq. (^) with respect 
to the disorder amplitude r. 

The maximum of the central charge is found in Fig. ^ 
at the amplitude of disorder r* ~ 3.881. The scaling 
dimension of the average order parameter, obtained at 
r*, is given in table ^ and compared to other disorder 
amplitudes. Again, there is a convincing agreement with 
the results obtained with the binary distribution. 
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With the continuous distribution at r = e, Olson and 
Young measured shghtly higher values for different val- 
ues of q. We propose here a possible explanation of 
the small discrepancy: It was shown that the larger 
the number of state q of the Potts model the stronger 
the disorder should be to reach the disordered fixed 
point regime ^4 19,5^]. Thus, for g > 3 random bond 
Potts models, the ideal disorder amplitude should be- 
come larger and larger, far from the value r = e ~ 2.718 
implicitly used by Olson and Young. On the other hand, 
since in the weak disorder regime the average exponent 
is continuously growing with the disorder amplitude r, a 
too weak randomness cannot be the explanation of the 
slightly too high values for the exponents obtained by 
Olson and Young. A possible origin of the small dis- 
agreement can be found in the ensemble average: If the 
number of disorder realizations is too small, the average 
behaviour will give an exponent closer to the typical one, 
and thus too large. 

TABLE V. Scaling dimension x'^ of the average order pa- 
rameter for the q — 3 random bond Potts model with a contin- 
uous distribution at the optimal disorder amplitude (written 
in bold face) and other amplitudes as well, compared to the 
results obtained in the binary case and for the pure model. 



InZ = lim„_ 



by introduction of n identical copies 



Distribution 


disorder amplitude 


<(3) 


Pure 


r = 1 


0.1333 


Binary 


r* ~ 5.363 


0.1347(11) 


Continuous 


r = 1.133 


0.1338(1) 




r = e 


0.1340(9) 




r* ~ 3.881 


0.1344(13) 




r = 7.389 


0.1348(17) 



B. Replica Symmetry 

The question of a possible breaking of replica symme- 
try |]57|| in disordered systems is very controversial and 
far from being settled, especially in spin glasses (see e.g 
Refs. p8|-|6l|). In the context of disordered Potts fer- 
romagnets, the question was first asked by Dotsenko et 
al ||9|. The Hamihonian (|l|) is rewritten 

- I3n[a] = Ko J2 + II (-^y ^ Ko)Sa,,a, 

(iJ) (id) 

= -pno[<j]-pn'[a] (31) 

where the deviation from the pure system can be written 
in the continuum limit 



(3n' 



(32) 



with t{x) — K{x) — Kc- The average free energy 
F = —ksThiZ can be obtained using the identity 



(labelled by an index a) of the model, coupled by their 
energy densities. After integration over a Gaussian prob- 
ability distribution centered on the value tq and with 
variance A, one is led to: 



Z" = Tr exp 



(a) 



To 



d^a; 



Eea(a 



ea{x)eb{x) 



(33) 



The first term governed by the average coupling tq pro- 
duces a shift in the critical temperature, while the sec- 
ond term couples the replicas with each other. In a 
replica symmetric scenario, the couplings A between 
the different copies of the model are identical, while in 
a Replica Symmetry Breaking scheme, these couplings 
are Parisi matrices and can take different values Aa?,. 
Treated as a perturbation this coupling term leads to 
different fixed point structures and finally to different 
scaling dimensions for the moments of the correlation 
functions. In order to test between Replica Symme- 
try and Replica Symmetry Breaking schemes, Dotsenko 
et al performed a second order expansion of the ex- 
ponent of the second moment of the spin-spin correla- 
tion function decay in both cases (Equations ( pd] ) and 
(|l^). Previous MC simulations have been performed at 
g = 3 but were not completely conclusive, although in 
favour of Replica Symmetry: The perturbation expan- 
sion leads to a;(.2(3) = 0.1176 and x^2(3) = 0.1201 ac- 

(|12| ), while previous numerical 
19r0.1140(5) ra, 0.116(1) 0] 



cording to Eq. (|ll| ) and 
results lead to 0.113(1) 
and 0.119(2) |l|. 

In this section, we report new conclusive results for 
different values of q. Close to g = 2, the proximity of 
the marginally irrelevant Ising FP will surely alter the 
data, as a reminiscent effect of the logarithmic correc- 
tions present exactly &i q = 2 for the second moment. 
Too large values of q on the other hand are not very 
helpful in order to check perturbation expansions which 
break down when one explores higher values of the ex- 
pansion parameter (as given for example in table ||) . One 
thus has to balance these two extreme situations and the 
comparison between numerical data and perturbation re- 
sults should be conclusive around q = i. The TM tech- 
nique thus appears to be well adapted, since it is capable 
to deal with non integer values of q. 

The comparison is shown in Fig. ^ for the bimodal 
probability distribution and the results are also given in 
table VI. In the convenient domain for the test, around 
5 = 3, results are written in bold face. The agreement 
with Replica Symmetry is quite convincing for the bi- 
modal probability distribution. We note that with the 
continuous distribution of Eq. ( p9| ) , we obtain also a very 
good value at g = 3: 0.1173(14) and with the ternary 
distribution at r = 8 we get 0.1182(12). 
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0.14 



0.13 



0.12 



0.11 



Pure system 




• • TM, bimodal distribution 


- 


RS 




RSB 


- 









2.5 



3.5 



q 



FIG. 7. Exponent of the second moment of the spin-spin 
correlation function as a function of the number of states of 
the disordered Potts model (binary disorder). The compar- 
ison is done with Reglica Symmetry and Replica Symmetry 
Breaking scenarios |49]. The agreement with the RS result 
is quite good around g = 3. When q is close to 2, the dis- 
crepancy can be attributed to the weak relevance of disorder. 
We indeed used a simple exponential fit as can be expected 
at a stable disordered FP, but at g = 2, one knows from 
Ludwig's results that logarithmic corrections must be added. 
These corrections can also influence the vicinity of g = 2 in a 
numerical approach. 



TABLE VI. Decay exponent of the second moment of the 
spin-spin correlation function compared to Replica Symmetry 
and Replica Symmetry Breaking expressions of Eqs. (^l|) and 
. The results written in bold face correspond to the range 
of values of g where the agreement is particularly satisfactory. 





Perturbative results 


TM result 


q 






Binary Disorder 


2.25 


0.12213 


0.12229 


0.1204(5) 


2.5 


0.12002 


0.12067 


0.1194(8) 


2.75 


0.11854 


0.11997 


0.1185(10) 


3. 


0.11761 


0.12011 


0.1177(12) 


3.25 


0.11718 


0.12110 


0.1172(14) 


3.5 


0.11723 


0.12304 


0.1169(16) 








Ternary Disorder 


3. 


0.11761 


0.12011 


0.1182(12) 








Continuous Disorder 


3. 


0.11761 


0.12011 


0.1173(14) 



C. Multifractality 



0.15 



0.1 



0.05 




• TM results 

1st order 

2nd order 



0.15 



a. 



0.1 



0.05 



q=?). 



0.15 



0.1 



0.05 



^=2.75 



'Tliv 



2 3 
P 



FIG. 8. Comparison of the multifractal exponents (reduced 

i/p 

moment of the correlation function {(j{Q)o{u))^ ) with the 
second order expansion of Lewis in the RS scheme p7| for dif- 
ferent values of q indicated in the figure (bimodal probability 
distribution) . 



The multiscaling behaviour of the spin-spin correlation 

functions is noticeable in the p— dependent set of expo- 

i/p 

nents of the reduced moments {a{0)a{p))P . The sec- 
ond moment has already been computed when looking 
for Replica Symmetry, and it can be generalized in the 
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strip geometry using Eq. (|T]) . We performed an exhaus- 
tive computation of 50 different moments in the range 
< p < 5 in the strip geometry, and the ass ociated scal- 
ing dimensions followed from a semi-log fit In {a{0)a{u))P 
vs Inw, according to Eq. (pT|), followed by an extrapola- 
tion to L — > oo. Examples for q = 2.75, 3 and 3.25 are 
shown in Fig. |^ (bimodal probability distribution) where 
the numerical results are also compared to the first order 
expansion of Ludwig and to the second order expansion 
in the RS scheme in Eq. (^0|) . The second order result is 
clearly very good up to values of p close to 3 and then 
breaks down as already noticed by Lewis |48) . 

An alternate presentation of the results (used e.g. by 
Ludwig ||] ) is given by the scaling dimension of the mo- 
ment of the correlation function itself, {a{0)a{p))P (not 

i/p 

the reduced function ((7(0)(7(/9))P ). The scaling dimen- 
sion is thus simply px'^p{q), hereafter denoted by X'^p(q) 
by a simple extension of Ludwig's notation. An exam- 
ple, with q = 3., is shown in Fig. ^ where we have also 
shown the results obtained with the continuous probabil- 
ity distribution at the optimal disorder amplitude. Once 
again, we find a promising agreement between the nu- 
merical data and the perturbative result which confirms 
universality. 



0.75 



q=3. 



0.5 



1st order 

2nd order 

• Bimodal distribution 

* Continuous distribution 



0.25 - 








FIG. 9. Comparison of the multifractal exponents (moment 
of the correlation function ((t(0)(7(u))p) with the second order 
expansion of Lewis in the RS scheme |M for both the bimodal 
and the continuous probability distributions. 

What is interesting in this latter presentation is the 
link with other fields where multifractality is observed. 
One then usually introduces a universal function, the 
multiscaling function H{a), which is simply the Legen- 
dre transform of the set of independent scaling indexes 
X'^p(q). Setting dX'^p{q) — adp, this function is simply 
obtained by H{a) — X'^p{q) — ap. The geometrical in- 
terpretation of this Legendre transform follows from the 

relation 4^ = — n where a is defined by ^'^"a ^'^^ = «■ 

aa ^ ^ op 



The scaling dimension x'^p {q) is obtained on the plot of 
H{a) by the intercept of the tangent of slope —p with 
the abscissa axis. An example of multiscaling function 
H[a) deduced from the numerical data with the bimodal 
probability distribution is shown in Fig. ^ for q — i 
and the line of slope —2, leading to the exponent of the 
second moment (See table ^) is also shown. For compar- 
ison, the function _ff(fc)(a) deduced from the fc— th order 
(in yn) perturbative results of Ludwig (fc = 1) and Lewis 
(fc = 2) (pl]| ) are also plotted. Other values of q are shown 
in Fig. 



0.235 



TM, bimodal distribution 




0.1177 



FIG. 10. Universal function -ff(Q) for different q = ?>. The 
functions H^^i^^a) and Hf2)[oi) deduced from Eq. ( px| ) at first 
or second order, respectively, are shown for comparison. 




FIG. 11. Universal function H[a) for different values of q 
(bimodal probability distribution). 
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D. Correlation function probability distribution 

In Ref . 1^ , Ludwig presented a remarquable discussion 
of the spin-spin correlation function probability distribu- 
tion. He showed how all the relevant information on the 
large distance behaviour is encoded in the multifractal 
function H{a). In this section, we follow Ludwig's ar- 
guments and report a numerical study of the correlation 
function probability distribution in the cylinder geome- 
try. 

According to the results of the previous section, the 
moments of the spin-spin correlation function along the 
strip asymptotically behaves as follows: 



GP{u) EE {a{0)a{u))P 



Bpe 



(34) 



and are defined in terms of the probability distribution 
V[G{u)]: 



GP{u) = / dG{u)V[G{u)]GP{u). 
Jq 



(35) 



Following Ludwig, we introduce the variable Y{u) = 
-lnG(w) and write GP{u) = e'P^'-'^l Using the iden- 
tity P[G{u)]dG = V[Y{u)]dY and equations ^) and 
(p5|), one obtains 



(36) 



which leads to the expression of the probability distribu- 
tion by inverting the Laplace transform ((5 > 0): 



V[Yiu)] 



1 
2i7r 



(5+ioo 



dp Bpe 



(37) 



Fol- 



The amplitude Bp is weakly depending on p. 
lowing Ludwig, it can be rewritten as Bp = 

exp — (~ 2Tr u^l ) ' '^^^ forgetten, since it only 
introduces a small correction when 2'ku/ L — > oo. Let us 
define the function h{p) = X'^p — 2^i/L -P- the large 
distance limit 2ttu/L oo, the integral can be evaluated 
by the saddle-point approximation at the minimum po of 

Hp): 



d'p^'^' 



Po 



Y{u) 
27ru/L 



(38) 



Instead of Y{u), we define the scaled variable a — '^^"^ 



27ru/L ' 

and the saddle point value at po only depends on this 
variable h{po) — H{a), where H{a) is nothing but the 



multifractal function defined in the previous section. We 
thus obtain the probability distribution 



V[Y{u)]^exp 



2nu f Yju) \ 
L \2ttu/LJ 



or, using 7'[y(u)]dF = V{a)da, 

2'KU 



, , 2lTU 

Vya) ~ — — exp 
ij 



L 



-H(a) 



(39) 



(40) 



The multifractal function contains the essential informa- 
tion on the probability distribution. In order to check 
this expression, the value of H{a) at fixed a is extracted 
by fitting the probability distribution to the expression 



In 7^(0;) = const 



, 2ttu 2'KU , , , ^ 

In— —Hia). (41) 



It is shown in Fig. 12 where the probability distribution 
of the spin-spin correlation function was obtained after 
collecting the results over 96 000 disorder realizations in 
50 classes. The values of H{a) are slightly too large, 
compared to the results presented in the previous sec- 
tion. We can indeed observe in Fig. ^ a deviation from 
the linear behaviour which would be expected with these 
variables, and the shorter the distance m, the larger the 
deviation. 

This could be due to a correction to the leading be- 
haviour given by the saddle-point approximation If 
we expand the function h{p) close to po, h(p) ~ H(a) + 
^h"{po){p — pq)^, with /i"(po) > we obtain, instead of 
Eq. (p9|), the following result for the probability distri- 
bution V[Y{u)] H: 



nviu)] ^ ( ^ 



-1/2 



exp 



27TU^ f Yju) \ 

L \2ku/l) 



(42) 



and a correction appears in V{ol): 

In-P(a) In = const H{a). (43) 

2 L L 

This is shown in Fig. |l^ where a linear behaviour is now 
obtained in the whole ran ge o f values ofu/L. A linear fit 
in the coordinates of Fig. |l3| gives the value of the mul- 
tifractal function H{a) which can be comp ared to the 
results of Fig. |lo| obtained in section IV C| 
performed for all values of 0.034 < a 



This fit is 
< 0.15 for (7 = 3 
for the cases of the bimodal probability distribution and 
of the continuous distribution. The results are shown in 
Fig. |l|. 



■^We mention here that a possible correction to the saddle- 
point approximation has been suggested by Olson and Young 
in their study of higher values of q. 
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FIG. 12. Behaviour of the probabihty distribution 'P{a) as 
a function of the distance along the strip 2nu/L. The rescaled 
position along the strip is written u = 2nu/L. Three fixed 
values of a are shown, the opened and filled symbols respec- 
tively correspond to the strip widths L = 6 and 7. A good 
collapse of the data at both sizes is observed, but the be- 
haviour displays a deviation from linearity at small distances 
{q = 3, binary distribution). 



FIG. 14. Multifractal function H{a) as it is deduced from 
the fit of the probability distribution Via) to Eq. ( |43| ) ac- 
counting for the correction near the saddle-point approxima- 
tion. It is compared to the results of Fig. |l^ in solid line. The 
insert shows slightly too large values of H{a) when deduced 
from Eq. @. 

This latter figure shows that the correlation function 
probability distribution is entirely determined by the uni- 
versal multifractal function H{a). 



I 



-6 - 



a=0.075 
''■■■««e".=*^.,a=0.100 



-11 



«fcoo 

a=0.125 



-16 



2nu/L 



50 



FIG. 13. Same as Fig. ^ accounting for the correction in 
Eq. (Usb close to the saddle-point approximation. 



V. CONCLUSIONS 

In this paper, we have investigated the critical be- 
haviour of the moments of the spin-spin correlation func- 
tions of two-dimensional random bond Potts ferromag- 
nets using Transfer Matrix techniques and conformal 
methods. New features of our present work are the fol- 
lowing. 

i) As far as we know, universality of the critical be- 
haviour of the moments of the correlation func- 
tion was checked for the first time in such sys- 
tems^. This statement follows from the numerical 
evidence that the scaling dimension of the average 
spin-spin correlation function, as well as those of 
higher moments or typical behaviour do not de- 
pend on the details of the probability distribution, 
provided that the computations are performed at 
the disordered fixed point given by the maximum 
condition of the central charge. The exponent of 
the average correlation function is furthermore in 



^We mention here that a MC study in the first-order regime 
of the pure model {q — 5) was recently reported wgere ran- 
dom bond disorder and dilution were found to belong to the 
same universality class |63|. 
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very good agreement with theoretical results using 
perturbative conformal field theory available in the 
literature. 

ii) The question of a possible breaking of Replica Sym- 

metry is also considered and the numerical data 
strongly support the absence of Replica Symmetry 
Breaking. The problem is investigated through the 
comparison of the scaling dimension of the second 
moment of the correlation function, which is com- 
pared to perturbative results obtained within both 
schemes. Although previous numerical results (us- 
ing Monte Carlo simulations) which led to similar 
conclusions were already reported at g = 3, we be- 
lieve that our study is conclusive, since it extends 
the work to other non integer values of q. 

iii) The multiscaling behaviour of the spin-spin corre- 

lation function is investigated. The exponential 
decay of the moments of the correlation functions 
along very long strips is used to deduce numerically 
the corresponding critical exponents. These dimen- 
sions continuously depend on the moment order, 
as a consequence of multifractality. They are fur- 
thermore very weakly dependent of the probability 
distribution. At low moment order, the numerical 
results are furthermore in agreement with a per- 
turbative result obtained within the Replica Sym- 
metry scheme. When the moment order increases, 
a discrepancy is observed, resulting from the lack 
of validity of the perturbation expansion, and pos- 
sibly of a numerical determination which becomes 
less precise. The multifractal function, given by 
the Legendre transform of the set of independent 
scaling dimensions is also computed for different 
values of q. It is shown that this universal function 
completely determines the shape of the correlation 
function probability distribution. 

The main result of this paper is probably to show that 
universality in random systems has to be understood in 
the sense of a critical behaviour which does not depend 
on the choice of the probability distribution (this is only 
true up to some extent, since special distributions which 
do not obey the central limit theorem, like Levy flights, 
would certainly lead to different results). By critical be- 
haviour, here we mean the behaviour of all the moments 
of a physical quantity, entirely contained in the multifrac- 
tal function or the correlation function probability distri- 
bution but we want to stress that lack of self-averaging 
does not imply absence of universality. 

We also note that logarithmic corrections were recently 
reported in the disconnected energy-energy correlations 
function and that disorder was shown to induce 
non-vanishing cross-correlations between spin-spin and 
energy-energy moments | ]65| . The multiscaling of energy 
correlations has been studied very recently by J.L. Ja- 
cobsen [M. 
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